EarthML Menu

Seattle Lidar

Right click to download this notebook from GitHub.

Visualize Lidar Scattered Point Elevation Data

This notebook uses Datashader to visualize Lidar elevation data from the Puget Sound Lidar consortium, a source of Lidar data for the Puget Sound region of the state of Washington, U.S.A.

Lidar Elevation Data

Example X,Y,Z scattered point elevation data from the unpacked 7zip files (unpacked as .gnd files)

! head ../data/q47122d2101.gnd
X,Y,Z
1291149.60,181033.64,467.95
1291113.29,181032.53,460.24
1291065.38,181035.74,451.41
1291113.16,181037.32,455.51
1291116.68,181037.42,456.20
1291162.42,181038.90,467.81
1291111.90,181038.15,454.89
1291066.62,181036.73,451.41
1291019.10,181035.20,451.64

The Seattle area example below loads 25 .gnd elevation files like the one above. We'll download, cache and read the data using intake.

In [1]:
import intake

cat = intake.open_catalog('../catalog.yml')
list(cat)
Out[1]:
['landsat_5_small',
 'landsat_8_small',
 'landsat_5',
 'landsat_8',
 'google_landsat_band',
 'amazon_landsat_band',
 'fluxnet_daily',
 'fluxnet_metadata',
 'seattle_lidar']
In [2]:
lidar = cat.seattle_lidar()
df = lidar.to_dask()
df.head()









Out[2]:












X
Y
Z






0
1293274.13
181399.80
435.79




1
1293278.38
181393.05
434.58




2
1293272.96
181401.08
435.86




3
1293271.02
181401.13
435.73




4
1293270.85
181407.80
436.15


























Geographic Metadata
Since the data are geographic, we need to know the coordinate reference system (CRS) of the X and Y. All the geographic metadata is stored in the data source entry in the intake catalog.












In [3]:






from pyproj import Proj, transform















In [4]:






lidar.metadata['crs']















Out[4]:




'State Plane Coordinate System Washington North FIPS 4601'
















In [5]:






washington_state_plane = Proj(init='epsg:2855')   # Washington State Plane North EPSG code
web_mercator = Proj(init='epsg:3857')             # Mercator projection EPSG code















In [6]:






FT_2_M = 0.3048  

def convert_coords(df):
    lon, lat = transform(washington_state_plane, web_mercator, df.X.values * FT_2_M, df.Y.values * FT_2_M)
    df['meterswest'], df['metersnorth'] = lon, lat 
    return df[['meterswest', 'metersnorth', 'Z']]



















Try out the convert_coords function on a subset of the data












In [7]:






convert_coords(df.head())















Out[7]:












meterswest
metersnorth
Z






0
-1.360741e+07
6.022196e+06
435.79




1
-1.360741e+07
6.022193e+06
434.58




2
-1.360741e+07
6.022197e+06
435.86




3
-1.360741e+07
6.022197e+06
435.73




4
-1.360741e+07
6.022200e+06
436.15


























Convert the coordinates
Since our real dataset is large and partitioned using dask, we need to think about how to apply the convert_coords function to our data.












In [8]:






import dask.distributed
import dask.delayed















In [9]:






dask.distributed.Client()















Out[9]:









Client





Cluster


    

  • Workers: 8
    • 

  • Cores: 8
    • 

  • Memory: 17.18 GB
    • 


























Explore the task graph to figure out a performant way to split up the coords conversion. First we'll try with using dask.delayed.












In [10]:






dask.delayed(convert_coords)(df) \
    .visualize()















Out[10]:























We can see that even though we thought dask.delayed would help, in actuality we would be requiring all the processes to be done first and then the conversion would happen on the whole dataset in one go. Another approach would be to use dask.map_partitions to do the conversion on each piece of the data.












In [11]:






df_merc = df.map_partitions(convert_coords)
df_merc.visualize()















Out[11]:























Now that we have set up the task graph, we can use df_merc directly to do the computations on the fly. However, since this dataset fits in memory on my computer, I will do the computation and keep the output in memory for quick use when plotting.


NOTE: This next cell takes about a minute to run. Take a look at the dask dashboard at the location specified above to watch the task progression.












In [12]:






%time dataset = df_merc.compute()



















CPU times: user 15.4 s, sys: 13 s, total: 28.3 s
Wall time: 1min 1s





















If all the data doesn't fit in memory on your machine, try downsampling the data from each file to only keep 1/100th of the total data. To avoid unnecessary computation, it is better to do the downsampling first and then convert the coordinates.












In [13]:






small = df.map_partitions(lambda x: convert_coords(x.sample(frac=0.01)))
small.visualize()















Out[13]:



















In [14]:






%time small_dataset = small.compute()



















CPU times: user 1.68 s, sys: 325 ms, total: 2.01 s
Wall time: 16.8 s





















Plot












In [15]:






import geoviews as gv
import holoviews as hv
from holoviews import opts
from holoviews.operation.datashader import datashade, rasterize, rd

hv.extension('bokeh')







































































To simplify the exploration of the time required to display different data, define a function that accepts a regular pandas dataframe or a dask delayed dataframe.












In [16]:






def plot(data, **kwargs):
    """Plot point elevation data, rasterizing by mean elevation"""
    points = hv.Points(data, kdims=['meterswest', 'metersnorth'], vdims=['Z'])
    image_opts = opts.Image(tools=['hover'], cmap='blues_r', colorbar=True,
                            width=800, height=800, xaxis=None, yaxis=None)
    return rasterize(points, aggregator=rd.mean('Z'), precompute=True, **kwargs).options(image_opts)



















We'll also declare a tiler to use for background imagery.












In [17]:






tiles = gv.tile_sources.EsriImagery()



















Then we will construct the plot out of rasterized point data and tiles.












In [18]:






tiles * plot(small_dataset)















Out[18]:

























Benchmarking












In [19]:






%time display(tiles * plot(small_dataset))
































CPU times: user 300 ms, sys: 33.7 ms, total: 334 ms
Wall time: 212 ms

















In [20]:






%time display(tiles * plot(dataset))
































CPU times: user 7.52 s, sys: 2.65 s, total: 10.2 s
Wall time: 7.27 s

















In [21]:






%time display(tiles * plot(df_merc))
































CPU times: user 9.86 s, sys: 1.27 s, total: 11.1 s
Wall time: 2min 19s





















Visualize Raster of Point Data
If we compute a raster of the point data then we don't need to store as much data in memory, which should allow faster interactions, and allow use with lots of other tools.












In [22]:






%time raster = plot(dataset, dynamic=False, width=1000, height=1000).data



















CPU times: user 6.56 s, sys: 2.02 s, total: 8.58 s
Wall time: 6.29 s





















The result is an xarray.Dataset with x and y coordinates and a 2D array of Z values.












In [23]:






raster.metersnorth[1] - raster.metersnorth[0]















Out[23]:




<xarray.DataArray 'metersnorth' ()>
array(10.292236)




















With these data we can use the geo tools in datashader to compute and visualize the elevation using hillshading for instance. See Datashader User Guide #8 for more datashader geo tooling.












In [24]:






from datashader.geo import hillshade















In [25]:






illuminated = hillshade(raster.Z)
hv_shaded = hv.Image((raster.meterswest, raster.metersnorth, illuminated))















In [26]:






tiles * rasterize(hv_shaded.options(opts.Image(cmap='blues', height=600, width=600)))















Out[26]:

























Ideas for future work
It'd be nice to have a rasterio writer from xarray so that we could easily write chunked geotiffs from xarray objects.


Something like:


raster.to_rasterio(path=None, mode='w', compute=True)